第9回 Pythonによる回帰分析の実装
y_i = \beta_0 + \beta_1 x_i + u_i
をモデルとして、回帰分析をする。
wage_and_productivity.xlsx| 2020年 | 労働生産性(円/時間) | 賃金所得(千円) | |
|---|---|---|---|
| 0 | 北海道 | 8997.702921 | 4232.5 |
| 1 | 青森 | 8946.026251 | 3667.9 |
| 2 | 岩手 | 8336.011457 | 3790.4 |
| 3 | 宮城 | 8875.070576 | 4459.4 |
| 4 | 秋田 | 8700.560140 | 3758.0 |
[[1.00000000e+00 8.99770292e+03]
[1.00000000e+00 8.94602625e+03]
[1.00000000e+00 8.33601146e+03]
[1.00000000e+00 8.87507058e+03]
[1.00000000e+00 8.70056014e+03]]
statsmodelsでOLSする場合は、sm.OLS()とfit()がセット。
| Dep. Variable: | y | R-squared: | 0.181 |
| Model: | OLS | Adj. R-squared: | 0.163 |
| Method: | Least Squares | F-statistic: | 9.978 |
| Date: | Fri, 12 Sep 2025 | Prob (F-statistic): | 0.00283 |
| Time: | 20:37:56 | Log-Likelihood: | -350.07 |
| No. Observations: | 47 | AIC: | 704.1 |
| Df Residuals: | 45 | BIC: | 707.8 |
| Df Model: | 1 | ||
| Covariance Type: | nonrobust |
| coef | std err | t | P>|t| | [0.025 | 0.975] | |
| const | 2681.0906 | 550.477 | 4.870 | 0.000 | 1572.373 | 3789.808 |
| x1 | 0.1830 | 0.058 | 3.159 | 0.003 | 0.066 | 0.300 |
| Omnibus: | 13.961 | Durbin-Watson: | 1.148 |
| Prob(Omnibus): | 0.001 | Jarque-Bera (JB): | 15.953 |
| Skew: | 1.089 | Prob(JB): | 0.000343 |
| Kurtosis: | 4.844 | Cond. No. | 8.45e+04 |
| 対応する値 | 回帰係数 | 標準誤差 | t統計量 | p値 | |
|---|---|---|---|---|---|
| 表示 | coef | std err | t | P>|t| | |
| const | 2681.1 | 550.5 | 4.870 | 0.000 | |
| x1 | 0.183 | 0.058 | 3.159 | 0.003 |
## statsmodelで回帰
## OLS()とfit()をセットで結果を出す
model = sm.OLS(y,X) # yをXにOLSで回帰
result = model.fit() # 先のmodelに.fit()をつけると、結果の計算がされる
#result = sm.OLS(y,X).fit() # こっちでもOK
beta = result.params # 回帰係数
SE = result.bse # 標準誤差
t = result.tvalues # t統計量
print("回帰係数は"+str(beta)+"です。")
print()
print("標準誤差は"+str(SE)+"です。")
print()
print("t統計量は"+str(t)+"です。")回帰係数は[2.68109056e+03 1.82966021e-01]です。
標準誤差は[5.50477090e+02 5.79230982e-02]です。
t統計量は[4.87048528 3.15877477]です。
デフォルトの回帰分析では、
であることを仮定して標準誤差を計算しています。
一方で、近年の分析では
fit()のオプションにcov_typeがあり、これを"HC0"と指定すれば不均一分散にも対応できるような標準誤差が計算できます。
## statsmodelで回帰
model = sm.OLS(y,X)
result0 = model.fit() # 均一分散を仮定
result1 = model.fit(cov_type="HC0") # 不均一分散に頑健
SE0 = result0.bse # 均一分散仮定のSE
SE1 = result1.bse # 不均一分散に頑健なSE
t0 = result0.tvalues
t1 = result1.tvalues
print("---均一分散を仮定すると---")
print("標準誤差は"+str(SE0))
print("t統計量は"+str(t0))
print()
print("---不均一分散に頑健に計算すると---")
print("標準誤差は"+str(SE1))
print("t統計量は"+str(t1))---均一分散を仮定すると---
標準誤差は[5.50477090e+02 5.79230982e-02]
t統計量は[4.87048528 3.15877477]
---不均一分散に頑健に計算すると---
標準誤差は[5.20905463e+02 5.52152006e-02]
t統計量は[5.14698108 3.31368934]
こういった推計をする際に、説明変数と被説明変数についてそれぞれ対数を取ってから回帰を行うことがあります。
\log Y = \beta_0 + \beta_1 \log X + u
## 変数を設定
x = np.log(df["労働生産性(円/時間)"].values) # 労働生産性の列の値をxとする
#x = np.log(df.iloc[:,1].values) # こっちでもOK
X = sm.add_constant(x) # 回帰式に定数項部分を足す
y = np.log(df["賃金所得(千円)"].values) # 賃金所得をyにする
print(X[:5]) # Xの最初の方を表示する
## OLS()とfit()をセットで結果を出す
model = sm.OLS(y,X) # yをXにOLSで回帰
result = model.fit(cov_type="HC0") # 先のmodelに.fit()をつけると、結果の計算がされる
# cov_type="HC0"を指定
#result = sm.OLS(y,X).fit(cov_type="HC0") # こっちでもOK
result.summary() # 結果を表示する。[[1. 9.10472459]
[1. 9.09896472]
[1. 9.02834014]
[1. 9.09100157]
[1. 9.07114269]]
| Dep. Variable: | y | R-squared: | 0.210 |
| Model: | OLS | Adj. R-squared: | 0.193 |
| Method: | Least Squares | F-statistic: | 13.41 |
| Date: | Fri, 12 Sep 2025 | Prob (F-statistic): | 0.000656 |
| Time: | 20:37:56 | Log-Likelihood: | 46.323 |
| No. Observations: | 47 | AIC: | -88.65 |
| Df Residuals: | 45 | BIC: | -84.95 |
| Df Model: | 1 | ||
| Covariance Type: | HC0 |
| coef | std err | z | P>|z| | [0.025 | 0.975] | |
| const | 4.5367 | 1.051 | 4.315 | 0.000 | 2.476 | 6.597 |
| x1 | 0.4208 | 0.115 | 3.662 | 0.000 | 0.196 | 0.646 |
| Omnibus: | 7.225 | Durbin-Watson: | 1.200 |
| Prob(Omnibus): | 0.027 | Jarque-Bera (JB): | 6.106 |
| Skew: | 0.774 | Prob(JB): | 0.0472 |
| Kurtosis: | 3.849 | Cond. No. | 765. |
回帰分析の説明変数が複数ある場合、この回帰分析を重回帰分析という。
この\beta_iについても、前回までの内容と同様に係数の推計と仮説検定が可能。
wage_and_others.xlsxをDL、pythonで読み込みをしてください。| 賃金所得(千円) | 労働生産性(円/時間) | 人口密度 | 生産年齢人口比率 | |
|---|---|---|---|---|
| 2020年 | ||||
| 北海道 | 4232.5 | 8997.702921 | 66.6 | 57.20614 |
| 青森 | 3667.9 | 8946.026251 | 128.3 | 55.72851 |
| 岩手 | 3790.4 | 8336.011457 | 79.2 | 55.41224 |
| 宮城 | 4459.4 | 8875.070576 | 316.1 | 60.18364 |
| 秋田 | 3758.0 | 8700.560140 | 82.4 | 52.83574 |
\begin{aligned} \log (賃金) = \beta_0 &+ \beta_1 \log(労働生産性) \\ &+ \beta_2 \log(人口密度) \\ &+ \beta_3 \log(生産年齢人口比率) + u \end{aligned}
で、帰無仮説H_0と対立仮説H_1をそれぞれ以下のように設定する。 H_0 :~ \beta_1 = 0 H_1 :~ \beta_1 \neq 0
ndarrayではなくpandasのDataFrameを指定することもできる。Xのチェック、最初の3行
const 労働生産性(円/時間) 人口密度 生産年齢人口比率
2020年
北海道 1.0 9.104725 4.198705 4.046661
青森 1.0 9.098965 4.854371 4.020492
岩手 1.0 9.028340 4.371976 4.014801
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.791
Model: OLS Adj. R-squared: 0.776
Method: Least Squares F-statistic: 121.2
Date: Fri, 12 Sep 2025 Prob (F-statistic): 5.34e-21
Time: 20:37:56 Log-Likelihood: 77.558
No. Observations: 47 AIC: -147.1
Df Residuals: 43 BIC: -139.7
Df Model: 3
Covariance Type: HC0
===============================================================================
coef std err z P>|z| [0.025 0.975]
-------------------------------------------------------------------------------
const 1.4936 0.766 1.949 0.051 -0.009 2.996
労働生産性(円/時間) 0.3491 0.078 4.494 0.000 0.197 0.501
人口密度 0.0406 0.011 3.696 0.000 0.019 0.062
生産年齢人口比率 0.8563 0.305 2.808 0.005 0.259 1.454
==============================================================================
Omnibus: 23.205 Durbin-Watson: 1.322
Prob(Omnibus): 0.000 Jarque-Bera (JB): 44.407
Skew: -1.393 Prob(JB): 2.28e-10
Kurtosis: 6.862 Cond. No. 1.94e+03
==============================================================================
Notes:
[1] Standard Errors are heteroscedasticity robust (HC0)
[2] The condition number is large, 1.94e+03. This might indicate that there are
strong multicollinearity or other numerical problems.
| 変数 | 係数 | 標準誤差 | z値(t値) | p値 |
|---|---|---|---|---|
| 表示 | coef | std err | z | P>|z| |
| const | 1.4936 | 0.766 | 1.949 | 0.051 |
| 労働生産性(円/時間) | 0.349 | 0.078 | 4.494 | 0.000 |
| 人口密度 | 0.0406 | 0.011 | 3.696 | 0.000 |
| 生産年齢人口比率 | 0.8563 | 0.305 | 2.808 | 0.005 |
DataFrameで説明変数を入れれば、変数にラベルを入れてくれる。第9回 Pythonによる回帰分析の実装